## "Public Preferences for International Law Compliance: 
##  Respecting Legal Obligations or Conforming to Common Practices?"
##
##
## Reproduce: 
##      - Figure 3
##
## Input:  
##      - data/data_cleaned.csv
##
## Output: 
##      - f3.pdf
##
## Saki Kuzushima
## 01/06/2023

source("code/functions_f3-4.R")
data <- read.csv(file="../data/data_cleaned.csv")

# output table
item <- c("surname3", "whale1", "hate1", "death1")
treat_group <- c("intl", "const", "oecd")
m_idx<- c(1, 0)
res <- expand.grid(m_idx=m_idx, item=item, treat=treat_group,
                   stringsAsFactors=F)
res$lower <- 0
res$point <- 0
res$upper <- 0
res$p <- 0
res$se <- 0
# moderator variable
data$id_cosmp_binary <- ifelse(data$id_cosmp <=3, 0, 1)
# NOTE: 
# 1 if people have stronger cosmopolitan identity. 0 otherise.

for (i in seq(1, nrow(res))){
  treatment <- paste0("treat_", res[i,"treat"])
  item <- res[i, "item"]
  m_idx <- res[i, "m_idx"]
  res[i, c("lower", "point", "upper", "p", "se")] <- 
    diffmean_block(data=data, treatment=treatment, item=item, 
             m="id_cosmp_binary", m_idx=m_idx)
}

# Multiple testing correction
res$p_bh <- p.adjust(res$p, method="BH")

# note: Based on BH, whaling and surname with constitutions, on the HI group  are significant 


# adaptive shrinkage
ash_fit <- ash(betahat=res$point, 
               sebetahat=res$se, 
               prior="uniform", 
               mixcompdist="unif",
               grange=c(-5, 5),
               outputlevel=2)
ash_result <- ash_fit$result
ash_ci <- ashci(ash_fit)
res_ash <- data.frame(m_idx=res$m_idx,
                      item=res$item,
                      treat=res$treat,
                      lower=ash_ci[,1], 
                      point=ash_result[,"PosteriorMean"],
                      upper=ash_ci[,2])
res_ash %>% filter(lower > 0)
# note: Based on Ash, whaling and surname with constitutions, on the HI group is signifncant
x_index_temp <- c(0.625, 0.875, 1.5, 1.5+0.625, 1.5+0.875, 3) 
x_index <- c(0, x_index_temp, 3+x_index_temp, 6+x_index_temp, 9+x_index_temp)


point_intl <- c(NA,res$point[1], res_ash$point[1], NA, 
                res$point[2], res_ash$point[2], NA, 
                res$point[3], res_ash$point[3], NA, 
                res$point[4], res_ash$point[4], NA,
                res$point[5], res_ash$point[5], NA, 
                res$point[6], res_ash$point[6], NA, 
                res$point[7], res_ash$point[7], NA,
                res$point[8], res_ash$point[8], NA)
upper_intl <- c(NA, res$upper[1], res_ash$upper[1], NA, 
                res$upper[2], res_ash$upper[2], NA, 
                res$upper[3], res_ash$upper[3], NA, 
                res$upper[4], res_ash$upper[4], NA,
                res$upper[5], res_ash$upper[5], NA, 
                res$upper[6], res_ash$upper[6], NA, 
                res$upper[7], res_ash$upper[7], NA,
                res$upper[8], res_ash$upper[8], NA)
lower_intl <- c(NA, res$lower[1], res_ash$lower[1], NA, 
                res$lower[2], res_ash$lower[2], NA, 
                res$lower[3], res_ash$lower[3], NA, 
                res$lower[4], res_ash$lower[4], NA,
                res$lower[5], res_ash$lower[5], NA, 
                res$lower[6], res_ash$lower[6], NA, 
                res$lower[7], res_ash$lower[7], NA,
                res$lower[8], res_ash$lower[8], NA)

point_const <- c(NA,res$point[9], res_ash$point[9], NA, 
                res$point[10], res_ash$point[10], NA, 
                res$point[11], res_ash$point[11], NA, 
                res$point[12], res_ash$point[12], NA,
                res$point[13], res_ash$point[13], NA, 
                res$point[14], res_ash$point[14], NA, 
                res$point[15], res_ash$point[15], NA,
                res$point[16], res_ash$point[16], NA)
upper_const <- c(NA,res$upper[9], res_ash$upper[9], NA, 
                res$upper[10], res_ash$upper[10], NA, 
                res$upper[11], res_ash$upper[11], NA, 
                res$upper[12], res_ash$upper[12], NA,
                res$upper[13], res_ash$upper[13], NA, 
                res$upper[14], res_ash$upper[14], NA, 
                res$upper[15], res_ash$upper[15], NA,
                res$upper[16], res_ash$upper[16], NA)
lower_const <- c(NA, res$lower[9], res_ash$lower[9], NA, 
                res$lower[10], res_ash$lower[10], NA, 
                res$lower[11], res_ash$lower[11], NA, 
                res$lower[12], res_ash$lower[12], NA,
                res$lower[13], res_ash$lower[13], NA, 
                res$lower[14], res_ash$lower[14], NA, 
                res$lower[15], res_ash$lower[15], NA,
                res$lower[16], res_ash$lower[16], NA)

point_oecd <- c(NA, res$point[17], res_ash$point[17], NA, 
                res$point[18], res_ash$point[18], NA, 
                res$point[19], res_ash$point[19], NA, 
                res$point[20], res_ash$point[20], NA,
                res$point[21], res_ash$point[21], NA, 
                res$point[22], res_ash$point[22], NA, 
                res$point[23], res_ash$point[23], NA,
                res$point[24], res_ash$point[24], NA)
upper_oecd <- c(NA, res$upper[17], res_ash$upper[17], NA, 
                res$upper[18], res_ash$upper[18], NA, 
                res$upper[19], res_ash$upper[19], NA, 
                res$upper[20], res_ash$upper[20], NA,
                res$upper[21], res_ash$upper[21], NA, 
                res$upper[22], res_ash$upper[22], NA, 
                res$upper[23], res_ash$upper[23], NA,
                res$upper[24], res_ash$upper[24], NA)
lower_oecd <- c(NA, res$lower[17], res_ash$lower[17], NA, 
                res$lower[18], res_ash$lower[18], NA, 
                res$lower[19], res_ash$lower[19], NA, 
                res$lower[20], res_ash$lower[20], NA,
                res$lower[21], res_ash$lower[21], NA, 
                res$lower[22], res_ash$lower[22], NA, 
                res$lower[23], res_ash$lower[23], NA,
                res$lower[24], res_ash$lower[24], NA)

point_style <- c(NA,1,16,NA,1,16,NA,2,17,NA,2,17,NA,0,15,NA,0,15,NA,5,18,NA,5,18,NA)
pdf("results/main/f3.pdf")
par(mfrow=c(3,1), mar=c(1,4,1,1)+0.1, oma=c(2,1,3,1)+0.1)
plot(x=x_index, y=point_intl, type="n",
     xlim=c(0.5, 11.5), ylim=c(-0.5, 0.5),
     pch=point_style, cex=1.5,
      yaxt="n", xaxt="n", xlab="", ylab="International Law", cex.lab=1.25,
     col="gray50")
abline(h=0, col="red")
abline(v=3, lty="dotted")
abline(v=6, lty="dotted")
abline(v=9, lty="dotted")
#abline(v=1.5, lty="dotted")
#abline(v=1.5+3, lty="dotted")
#abline(v=1.5+6, lty="dotted")
#abline(v=1.5+9, lty="dotted")
axis(side=3, at=c(1.5, 1.5+3, 1.5+6, 1.5+9),
     labels=c("Surname", "Whale", "Hate Speech", "Death Penalty"),
     tick=F,  line=0, outer=T, cex.axis=1.25)
axis(side=3, at=c(0.75, 2.25, 0.75+3, 2.25+3, 0.75+6, 2.25+6, 0.75+9, 2.25+9),
     labels=rep(c("HI", "LO"),4),
     tick=F,  line=0, cex.axis=1.15)
axis(side=2, at=seq(-0.5, 0.5, 0.25), labels=seq(-0.5, 0.5, 0.25), tick=T)
points(x=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       y=point_intl[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       pch=point_style[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       cex=c(rep(1.5,12), 1.5,2,1.5,2),
       col="grey50") 
segments(x0=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)],
         x1=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         y0=lower_intl[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         y1=upper_intl[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         cex=1.5,
         col="grey50") 
points(x=x_index[14],
       y=point_intl[14],
       pch=point_style[14],
       cex=1.5,
       col="black")
segments(x0=x_index[14],
         x1=x_index[14],
         y0=lower_intl[14],
         y1=upper_intl[14],
         cex=1.5,
         col="black")
plot(x=x_index, y=point_const, type="n",
     xlim=c(0.5, 11.5), ylim=c(-0.5, 0.5),
     pch=point_style, cex=1.5,
      yaxt="n", xaxt="n", xlab="", ylab="Constitution", cex.lab=1.25,
     col="gray50")
abline(h=0, col="red")
abline(v=3, lty="dotted")
abline(v=6, lty="dotted")
abline(v=9, lty="dotted")
#abline(v=1.5, lty="dotted")
#abline(v=1.5+3, lty="dotted")
#abline(v=1.5+6, lty="dotted")
#abline(v=1.5+9, lty="dotted")
axis(side=2, at=seq(-0.5, 0.5, 0.25), labels=seq(-0.5, 0.5, 0.25), tick=T)
points(x=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       y=point_const[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       pch=point_style[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       cex=c(rep(1.5,12), 1.5,2,1.5,2),
       col="grey50") 
segments(x0=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)],
         x1=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         y0=lower_const[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         y1=upper_const[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         cex=1.5,
         col="grey50") 
points(x=x_index[c(2,3,8,9,14)],
       y=point_const[c(2,3,8,9,14)],
       pch=point_style[c(2,3,8,9,14)],
       cex=1.5,
       col="black")
segments(x0=x_index[c(2,3,8,9,14)],
         x1=x_index[c(2,3,8,9,14)],
         y0=lower_const[c(2,3,8,9,14)],
         y1=upper_const[c(2,3,8,9,14)],
         cex=1.5,
         col="black")
text(x=1.1, y=0.25, labels="BH"~symbol("\326"))
text(x=4.1, y=0.25, labels="BH"~symbol("\326"))
plot(x=x_index, y=point_oecd, type="n",
     xlim=c(0.5, 11.5), ylim=c(-0.5, 0.5),
     pch=point_style, cex=1.5,
      yaxt="n", xaxt="n", xlab="", ylab="Practices", cex.lab=1.25,
     col="gray50")
abline(h=0, col="red")
abline(v=3, lty="dotted")
abline(v=6, lty="dotted")
abline(v=9, lty="dotted")
#abline(v=1.5, lty="dotted")
#abline(v=1.5+3, lty="dotted")
#abline(v=1.5+6, lty="dotted")
#abline(v=1.5+9, lty="dotted")
axis(side=2, at=seq(-0.5, 0.5, 0.25), labels=seq(-0.5, 0.5, 0.25), tick=T)
points(x=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       y=point_oecd[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       pch=point_style[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
       cex=c(rep(1.5,12), 1.5,2,1.5,2),
       col="grey50") 
segments(x0=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)],
         x1=x_index[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         y0=lower_oecd[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         y1=upper_oecd[c(c(2,3,5,6),c(2,3,5,6)+6,c(2,3,5,6)+12,c(2,3,5,6)+18)], 
         cex=1.5,
         col="grey50") 
dev.off()
